##### R Code for 'Population Dynamics' #####

# https://compadre-db.org/Education/article/population-dynamics

# Compiled on 9 February 2021

##### Preliminaries #####

# Use the command install.packages("popdemo") if you haven't already
# downloaded the package

library(popdemo) # Activates the package once installed in your library

##### Deterministic models #####

# Add data
data(Tort)
Tort

# Choose a vector
vec <- runif(8)

# Project to 50 time intervals
Tortp <- project(Tort, vec, time = 50, standard.vec = TRUE)
Tortp

# Observe population vectors
vec(Tortp)[1:3,] # Here, we've limited the display to rows 1-3

# Plot population
plot(Tortp, log ="y")

##### Asymptotic dynamics #####

# Calculate dominant eigenvalue and eigenvectors
eigs(Tort, "all")

# Project to 50 time intervals with vector 'w'
Tortw <- eigs(Tort, "ss")
Tortpw <- project(Tort, Tortw, time = 50)

# Plot population
plot(Tortp, log = "y")
lines(0:50, Tortpw, lty = 2) # This line adds to the plot our recently projected population with the stable ratios

##### Transient dynamics #####

# Plot standardizations
Tortps <- project(Tort, vec, time = 50, 
                  standard.A = TRUE, standard.vec = TRUE)

Tortpws <- project(Tort, Tortw, time = 50, 
                   standard.A = TRUE, standard.vec = TRUE)

plot(Tortps, log = "y")
lines(Tortpws, lty = 2)

# Calculate damping ratio
dr(Tort)

##### Transient indices #####

# Create sample populations
TortAMP <- c(1, 1, 2, 3, 5, 8, 13, 21) # A population that amplifies
TortATT <- c(21, 13, 8, 5, 3, 2, 1, 1) # A population that attenuates
TortBTH <- c(0, 0, 0, 1, 0, 0, 0, 0) # A population that does both

# Bind them together
Tortvec3 <- cbind(AMP = TortAMP, ATT = TortATT, BTH = TortBTH)
Tortvec3

# Project to 50 time intervals
Tortp3 <- project(Tort, Tortvec3, time = 50, 
                  standard.A = TRUE, standard.vec = TRUE)

# Define points
max3 <- apply(Tortvec3[,c(1,3)], 2, maxamp, A = Tort)

max3t <- apply(Tortvec3[,c(1,3)], 2, function(x){
  maxamp(vector = x, A = Tort, return.t = TRUE)$t})

min3 <- apply(Tortvec3[,c(2,3)], 2, maxatt, A = Tort)

min3t <- apply(Tortvec3[,c(2,3)], 2, function(x){
  maxatt(vector = x, A = Tort, return.t = TRUE )$t})

# Plot population
# And finally, let's plot
plot(Tortp3, log = "y"); lines(Tortpws, lty = 2)
points(c(max3t, min3t),
       c(max3, min3),
       pch = 3, col = "red")